﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using AdaptiveDeNet.CSV;
using AdaptiveDeNet.RAND;
using System.Threading;
using System.Threading.Tasks;
using AdaptiveDeNet.DifferentialEvolution.Common;
using System.Diagnostics;

namespace AdaptiveDeNet.DifferentialEvolution.MultiThreaded
{
    public partial class Optimizer
    {
        //b6e6rl variables

        private const double H = 12.0;
        public long[] strategySuccesses;
        public double[] stratPropabilities = new double[12];
        private double n0;
        private double delta;

        /// <summary>
        /// b6e6rl competitive DE, tries to find global minimum point of given function
        /// </summary>
        /// <param name="function">objective function</param>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">problem dimesion</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="n0">controls changes in probability of strategy selection</param>
        /// <param name="delta">critial minimal value of strategy propability</param>
        /// <returns></returns>
        public CEC14Result b6e6rl(double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double n0, double delta, Queue<double> FESMultipliers, double fX, int funcNum)
        {
            CEC14Result result = new CEC14Result(); //retval

            mre = new ManualResetEvent(false);
            countdown = new CountdownEvent(threadCount);
            outerCountdown = new CountdownEvent(threadCount);

            lockers = new ManualResetEventSlim[threadCount];

            mreslis = new ManualResetEventSlim[threadCount];

            condition = true;

            this.n0 = n0;
            this.delta = delta;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i] = new ManualResetEventSlim(false);
            }


            for (int i = 0; i < mreslis.Length; i++)
            {
                mreslis[i] = new ManualResetEventSlim(false);
            }

            //random stream initialization
            RngStream stream = new RngStream();
            long seed = DateTime.Now.Ticks & Int32.MaxValue;
            RngStream.setPackageSeed(new long[] { seed, seed, seed, seed, seed, seed });

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            double[] uBlBDiff = new double[dim]; //upperBound - lowerBound 

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            int popEvals = N;

            int padding = 16;
            int spacing = 16;

            //pointers
            int[] populationPointers = new int[N];
            int[] fitnessPointers = new int[N];

            double[] mainPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] trialsPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];

            int helper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;

                double[] temp = new double[dim];
                double fEvalInner = 0;

                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    //member vector initialization
                    int start = padding + i * (dim + 1) + threadActual * spacing;
                    int end = start + dim;

                    populationPointers[helper] = start;
                    fitnessPointers[helper] = end;

                    for (int j = start, forDiff = 0; j < end; j++, forDiff++)
                    {
                        temp[forDiff] = nRand.NextDouble() * uBlBDiff[forDiff] + lowerBound[forDiff];
                    }

                    fEvalInner = fEval.Feval(temp); //evaluation
                    if (fEvalInner < Fmin) Fmin = fEvalInner; //is fMin ?

                    Array.Copy(temp, 0, mainPopHolder, start, dim); // temp -> mainPopHolder
                    mainPopHolder[end] = fEvalInner; // fEval -> mainPopHolder

                    helper++;
                }
            }

            Fmax = this.FindMax(mainPopHolder, fitnessPointers);
            Fmin = this.FindMin(mainPopHolder, fitnessPointers);

            // init strategy pool input variables, find real roots of polynomical equantion etc.
            // this should be a list of tuple variables F/Cr

            //get strategy F and Cr for binomical/exponantial crossover

            double[][] controlParameters = new double[12][];

            controlParameters = GetControlParameters(dim);

            //number strategy successes initialiazation/
            int longPadding = 16;
            int longSpacing = 16;

            strategySuccesses = new long[longPadding + (12 * longSpacing)];
            int[] stratPtrs = new int[12];

            for (int i = 0; i < 12; i++)
            {
                int offset = 16 + i * 16;

                strategySuccesses[offset] = 1;
                stratPtrs[i] = offset;
            }

            ResetStrategySelectionProb();
            ResetStrategySuccesses();

            long[] affinity = new long[threadCount];

            affinity[0] = 0x0001;
            affinity[1] = 0x0004;


            // main loop
            for (int k = 0; k < threadCount; k++)
            {
                int iThread = k;

                ThreadPool.QueueUserWorkItem((st) =>
                {
                    lockers[iThread] = new ManualResetEventSlim(true);


                    int start = N * iThread / threadCount;
                    int end = N * (iThread + 1) / threadCount;
                    int n1 = N - 1;
                    int dim1 = dim - 1;
                    RngStream localStream = new RngStream();

                    uint unmanagedId = GetCurrentThreadId();

                    ProcessThread myThread = (from ProcessThread entry in Process.GetCurrentProcess().Threads
                                              where entry.Id == unmanagedId
                                              select entry).First();

                    myThread.ProcessorAffinity = (IntPtr)affinity[iThread];



                    CEC14Lib.CEC14func localFunc = new CEC14Lib.CEC14func(30, funcNum);

                    double[] temp = new double[dim];
                    double innerFEval = 0;
                    int i = start;

                    //roulette
                    double rd = 0.0;
                    double sum = 0.0;

                    int strategyIndex = 0;
                    int iPtr = 0;

                    countdown.Signal();

                    do
                    {

                        lockers[iThread].Wait();
                        lockers[iThread].Reset();

                        for (i = start; i < end; i++)
                        {
                            //roulette
                            rd = localStream.NextDouble();
                            sum = 0.0;

                            strategyIndex = 0;
                            iPtr = 0;

                            do
                            {
                                sum += stratPropabilities[iPtr];
                                strategyIndex = iPtr;
                                iPtr++;

                            } while (sum < rd && iPtr < 12);

                            double F = controlParameters[strategyIndex][0];
                            double Cr = controlParameters[strategyIndex][1];

                            int r0Index = localStream.randInt(0, n1);
                            int r1Index = localStream.randInt(0, n1);
                            int r2Index = localStream.randInt(0, n1);

                            while (i == r0Index) r0Index = localStream.randInt(0, n1);

                            while (i == r1Index || r0Index == r1Index) r1Index = localStream.randInt(0, n1);

                            while (i == r2Index || r0Index == r2Index || r1Index == r2Index) r2Index = localStream.randInt(0, n1);

                            int jrand = localStream.randInt(0, dim1);

                            //RANDOM LOCATION
                            int pom = 0;

                            if (mainPopHolder[fitnessPointers[r2Index]] > mainPopHolder[fitnessPointers[r1Index]])
                            {
                                if (mainPopHolder[fitnessPointers[r1Index]] < mainPopHolder[fitnessPointers[r0Index]])
                                {
                                    pom = r0Index;
                                    r0Index = r1Index;
                                    r1Index = pom;
                                }
                            }
                            else
                            {
                                if (mainPopHolder[fitnessPointers[r2Index]] < mainPopHolder[fitnessPointers[r0Index]])
                                {
                                    pom = r0Index;
                                    r0Index = r2Index;
                                    r2Index = pom;
                                }
                            }

                            int innerStart = padding + i * (dim + 1) + iThread * spacing;
                            int innerEnd = innerStart + dim;

                            int popPtr = i;

                            if (strategyIndex <= 6)//binomical crossover
                            {
                                for (int j = innerStart, jPtr = 0; j < innerEnd; j++, jPtr++)
                                {
                                    //krizeni
                                    if (localStream.randU01() <= Cr || jPtr == jrand)
                                    {
                                        temp[jPtr] = mainPopHolder[populationPointers[r0Index] + jPtr] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jPtr] -
                                            mainPopHolder[populationPointers[r2Index] + jPtr]);

                                        if (temp[jPtr] < lowerBound[jPtr] || temp[jPtr] > upperBound[jPtr])
                                        {
                                            temp[jPtr] = Common.Common.Zrdcad(temp[jPtr],
                                                lowerBound[jPtr], upperBound[jPtr]);
                                        }
                                    }
                                    else
                                    {
                                        temp[jPtr] = mainPopHolder[j];
                                    }
                                }
                            }
                            else
                            {
                                int jRandtemp = jrand;

                                if (strategyIndex > 6)
                                {
                                    for (int o = 0; o < temp.Length; o++)
                                    {
                                        temp[o] = 0;
                                    }

                                    do
                                    {
                                        // 
                                        // x = r0 + F * (r1 - r2) - mutant
                                        temp[jRandtemp] = mainPopHolder[populationPointers[r0Index] + jRandtemp] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jRandtemp] -
                                            mainPopHolder[populationPointers[r2Index] + jRandtemp]);

                                        //zrdcadleni
                                        if (temp[jRandtemp] < lowerBound[jRandtemp] || temp[jRandtemp] > upperBound[jRandtemp])
                                        {
                                            //newVal = nRand.NextDouble() * (upperBound[jrandTemp] - lowerBound[jrandTemp]) + lowerBound[jrandTemp]; 
                                            temp[jRandtemp] = Common.Common.Zrdcad(temp[jRandtemp], lowerBound[jRandtemp], upperBound[jRandtemp]);
                                        }

                                        ++jRandtemp; //increment

                                    } while (localStream.NextDouble() <= Cr && jRandtemp < dim);

                                    //pokud rand <= Cr pretece populaci
                                    if (jRandtemp == dim - 1)
                                    {
                                        jRandtemp = 0;

                                        do
                                        {
                                            // 
                                            // x = r0 + F * (r1 - r2) - mutant
                                            temp[jRandtemp] = mainPopHolder[populationPointers[r0Index] + jRandtemp] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jRandtemp] -
                                            mainPopHolder[populationPointers[r2Index] + jRandtemp]);

                                            //zrdcadleni
                                            if (temp[jRandtemp] < lowerBound[jRandtemp] || temp[jRandtemp] > upperBound[jRandtemp])
                                            {
                                                //newVal = nRand.NextDouble() * (upperBound[jrandTemp] - lowerBound[jrandTemp]) + lowerBound[jrandTemp]; 
                                                temp[jRandtemp] = Common.Common.Zrdcad(temp[jRandtemp], lowerBound[jRandtemp], upperBound[jRandtemp]);
                                            }

                                            ++jRandtemp; //increment  

                                        } while (localStream.NextDouble() <= Cr && jRandtemp < jrand);

                                        if (jRandtemp < jrand)
                                        {
                                            for (int j = jRandtemp; j < jrand; j++)
                                            {
                                                temp[j] = mainPopHolder[innerStart + j];
                                            }
                                        }
                                    }
                                    else //jinak dokopirovani
                                    {
                                        for (int j = jRandtemp; j < dim; j++)
                                        {
                                            temp[j] = mainPopHolder[innerStart + j];
                                        }

                                        for (int j = 0; j < jrand; j++)
                                        {
                                            temp[j] = mainPopHolder[innerStart + j];
                                        }
                                    }//endElse
                                }//endIf
                            }//endElse



                            //Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);
                            //Array.Copy(trialsPopHolder, innerStart, temp, 0, dim);
                            //Buffer.BlockCopy(Q, innerStart, temp, 0, dim * sizeof(double));

                            innerFEval = localFunc.FEval(temp, 30, funcNum);
                            //Interlocked.Increment(ref truePopEvals);
                            //trialsFevals[i] = innerFEval;
                            //trialsPopHolder[innerEnd] = innerFEval;

                            if (innerFEval <= mainPopHolder[innerEnd])
                            {
                                Interlocked.Increment(ref strategySuccesses[stratPtrs[strategyIndex]]);

                                trialsPopHolder[innerEnd] = innerFEval;
                                Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);
                            }//succesful trial
                            else
                            {
                                trialsPopHolder[innerEnd] = mainPopHolder[innerEnd];
                                Array.Copy(mainPopHolder, innerStart, trialsPopHolder, innerStart, dim);
                            }

                        }


                        mreslis[iThread].Set();
                        //countdown.Signal();

                    } while (condition);

                    outerCountdown.Signal();

                });
            }

            countdown.Wait();

            while (popEvals < maxEvals)
            {


                for (int i = 0; i < threadCount; i++)
                {
                    lockers[i].Set();
                }

                for (int i = 0; i < mreslis.Length; i++)
                {

                    mreslis[i].Wait();
                    mreslis[i].Reset();
                }

                double[] pom;
                pom = mainPopHolder;
                mainPopHolder = trialsPopHolder;
                trialsPopHolder = pom;

                Fmax = this.FindMax(mainPopHolder, fitnessPointers);
                Fmin = this.FindMin(mainPopHolder, fitnessPointers);

                popEvals += N;

                EvalStrategyPool(stratPtrs);

                if (popEvals == (FESMultipliers.Peek() * maxEvals))
                {
                    result.Records.Add(Fmin - fX);
                    FESMultipliers.Dequeue();
                }

                //countdown.Reset(2);
            }

            condition = false;
            //Monitor.PulseAll(locker);
            //countdown.Wait();

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i].Set();
            }

            outerCountdown.Wait();

            result.FuncEvals = popEvals;


            return result;
        }

        /// <summary>
        /// b6e6rl competitive DE, tries to find global minimum point of given function
        /// </summary>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">problem dimesion</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="n0">controls changes in probability of strategy selection</param>
        /// <param name="delta">critial minimal value of strategy propability</param>
        /// <returns></returns>
        public Result b6e6rl(double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double n0, double delta)
        {
            mre = new ManualResetEvent(false);
            countdown = new CountdownEvent(threadCount);
            outerCountdown = new CountdownEvent(threadCount);

            lockers = new ManualResetEventSlim[threadCount];

            mreslis = new ManualResetEventSlim[threadCount];

            condition = true;

            this.n0 = n0;
            this.delta = delta;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i] = new ManualResetEventSlim(false);
            }


            for (int i = 0; i < mreslis.Length; i++)
            {
                mreslis[i] = new ManualResetEventSlim(false);
            }


            //random stream initialization
            RngStream stream = new RngStream();
            long seed = DateTime.Now.Ticks & Int32.MaxValue;
            RngStream.setPackageSeed(new long[] { seed, seed, seed, seed, seed, seed });

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            double[] uBlBDiff = new double[dim]; //upperBound - lowerBound 

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            int popEvals = N;

            int padding = 16;
            int spacing = 16;

            //pointers
            int[] populationPointers = new int[N];
            int[] fitnessPointers = new int[N];

            double[] mainPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] trialsPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];

            int helper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;

                double[] temp = new double[dim];
                double fEvalInner = 0;

                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    //member vector initialization
                    int start = padding + i * (dim + 1) + threadActual * spacing;
                    int end = start + dim;

                    populationPointers[helper] = start;
                    fitnessPointers[helper] = end;

                    for (int j = start, forDiff = 0; j < end; j++, forDiff++)
                    {
                        temp[forDiff] = nRand.NextDouble() * uBlBDiff[forDiff] + lowerBound[forDiff];
                    }

                    fEvalInner = fEval.Feval(temp); //evaluation
                    if (fEvalInner < Fmin) Fmin = fEvalInner; //is fMin ?

                    Array.Copy(temp, 0, mainPopHolder, start, dim); // temp -> mainPopHolder
                    mainPopHolder[end] = fEvalInner; // fEval -> mainPopHolder

                    helper++;
                }
            }

            Fmax = this.FindMax(mainPopHolder, fitnessPointers);
            Fmin = this.FindMin(mainPopHolder, fitnessPointers);

            // init strategy pool input variables, find real roots of polynomical equantion etc.
            // this should be a list of tuple variables F/Cr

            //get strategy F and Cr for binomical/exponantial crossover

            double[][] controlParameters = new double[12][];

            controlParameters = GetControlParameters(dim);

            //number strategy successes initialiazation/
            int longPadding = 16;
            int longSpacing = 16;

            strategySuccesses = new long[longPadding + (12 * longSpacing)];
            int[] stratPtrs = new int[12];

            for (int i = 0; i < 12; i++)
            {
                int offset = 16 + i * 16;

                strategySuccesses[offset] = 1;
                stratPtrs[i] = offset;
            }

            ResetStrategySelectionProb();
            ResetStrategySuccesses();

            // main loop
            for (int k = 0; k < threadCount; k++)
            {
                int iThread = k;

                ThreadPool.QueueUserWorkItem((st) =>
                {
                    lockers[iThread] = new ManualResetEventSlim(true);


                    int start = N * iThread / threadCount;
                    int end = N * (iThread + 1) / threadCount;
                    int n1 = N - 1;
                    int dim1 = dim - 1;
                    RngStream localStream = new RngStream();

                    double[] temp = new double[dim];
                    double innerFEval = 0;
                    int i = start;

                    //roulette
                    double rd = 0.0;
                    double sum = 0.0;

                    int strategyIndex = 0;
                    int iPtr = 0;

                    countdown.Signal();

                    do
                    {

                        lockers[iThread].Wait();
                        lockers[iThread].Reset();

                        for (i = start; i < end; i++)
                        {
                            //roulette
                            rd = localStream.NextDouble();
                            sum = 0.0;

                            strategyIndex = 0;
                            iPtr = 0;

                            do
                            {
                                sum += stratPropabilities[iPtr];
                                strategyIndex = iPtr;
                                iPtr++;

                            } while (sum < rd && iPtr < 12);

                            double F = controlParameters[strategyIndex][0];
                            double Cr = controlParameters[strategyIndex][1];

                            int r0Index = localStream.randInt(0, n1);
                            int r1Index = localStream.randInt(0, n1);
                            int r2Index = localStream.randInt(0, n1);

                            while (i == r0Index) r0Index = localStream.randInt(0, n1);

                            while (i == r1Index || r0Index == r1Index) r1Index = localStream.randInt(0, n1);

                            while (i == r2Index || r0Index == r2Index || r1Index == r2Index) r2Index = localStream.randInt(0, n1);

                            int jrand = localStream.randInt(0, dim1);

                            //RANDOM LOCALIZATION
                            int pom = 0;

                            if (mainPopHolder[fitnessPointers[r2Index]] > mainPopHolder[fitnessPointers[r1Index]])
                            {
                                if (mainPopHolder[fitnessPointers[r1Index]] < mainPopHolder[fitnessPointers[r0Index]])
                                {
                                    pom = r0Index;
                                    r0Index = r1Index;
                                    r1Index = pom;
                                }
                            }
                            else
                            {
                                if (mainPopHolder[fitnessPointers[r2Index]] < mainPopHolder[fitnessPointers[r0Index]])
                                {
                                    pom = r0Index;
                                    r0Index = r2Index;
                                    r2Index = pom;
                                }
                            }

                            int innerStart = padding + i * (dim + 1) + iThread * spacing;
                            int innerEnd = innerStart + dim;

                            int popPtr = i;

                            if (strategyIndex <= 6)//binomical crossover
                            {
                                for (int j = innerStart, jPtr = 0; j < innerEnd; j++, jPtr++)
                                {
                                    //krizeni
                                    if (localStream.randU01() <= Cr || jPtr == jrand)
                                    {
                                        temp[jPtr] = mainPopHolder[populationPointers[r0Index] + jPtr] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jPtr] -
                                            mainPopHolder[populationPointers[r2Index] + jPtr]);

                                        if (temp[jPtr] < lowerBound[jPtr] || temp[jPtr] > upperBound[jPtr])
                                        {
                                            temp[jPtr] = Common.Common.Zrdcad(temp[jPtr],
                                                lowerBound[jPtr], upperBound[jPtr]);
                                        }
                                    }
                                    else
                                    {
                                        temp[jPtr] = mainPopHolder[j];
                                    }
                                }
                            }
                            else
                            {
                                int jRandtemp = jrand;

                                if (strategyIndex > 6)
                                {
                                    for (int o = 0; o < temp.Length; o++)
                                    {
                                        temp[o] = 0;
                                    }

                                    do
                                    {
                                        // 
                                        // x = r0 + F * (r1 - r2) - mutant
                                        temp[jRandtemp] = mainPopHolder[populationPointers[r0Index] + jRandtemp] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jRandtemp] -
                                            mainPopHolder[populationPointers[r2Index] + jRandtemp]);

                                        //zrdcadleni
                                        if (temp[jRandtemp] < lowerBound[jRandtemp] || temp[jRandtemp] > upperBound[jRandtemp])
                                        {
                                            temp[jRandtemp] = Common.Common.Zrdcad(temp[jRandtemp], lowerBound[jRandtemp], upperBound[jRandtemp]);
                                        }

                                        ++jRandtemp; //increment

                                    } while (localStream.NextDouble() <= Cr && jRandtemp < dim);

                                    //pokud rand <= Cr pretece populaci
                                    if (jRandtemp == dim - 1)
                                    {
                                        jRandtemp = 0;

                                        do
                                        {
                                            // 
                                            // x = r0 + F * (r1 - r2) - mutant
                                            temp[jRandtemp] = mainPopHolder[populationPointers[r0Index] + jRandtemp] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jRandtemp] -
                                            mainPopHolder[populationPointers[r2Index] + jRandtemp]);

                                            //zrdcadleni
                                            if (temp[jRandtemp] < lowerBound[jRandtemp] || temp[jRandtemp] > upperBound[jRandtemp])
                                            {
                                                temp[jRandtemp] = Common.Common.Zrdcad(temp[jRandtemp], lowerBound[jRandtemp], upperBound[jRandtemp]);
                                            }

                                            ++jRandtemp; //increment  

                                        } while (localStream.NextDouble() <= Cr && jRandtemp < jrand);

                                        if (jRandtemp < jrand)
                                        {
                                            for (int j = jRandtemp; j < jrand; j++)
                                            {
                                                temp[j] = mainPopHolder[innerStart + j];
                                            }
                                        }
                                    }
                                    else //jinak dokopirovani
                                    {
                                        for (int j = jRandtemp; j < dim; j++)
                                        {
                                            temp[j] = mainPopHolder[innerStart + j];
                                        }

                                        for (int j = 0; j < jrand; j++)
                                        {
                                            temp[j] = mainPopHolder[innerStart + j];
                                        }
                                    }//endElse
                                }//endIf
                            }//endElse

                            innerFEval = fEval.Feval(temp);

                            if (innerFEval <= mainPopHolder[innerEnd])
                            {
                                Interlocked.Increment(ref strategySuccesses[stratPtrs[strategyIndex]]);

                                trialsPopHolder[innerEnd] = innerFEval;
                                Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);
                            }//succesful trial
                            else
                            {
                                trialsPopHolder[innerEnd] = mainPopHolder[innerEnd];
                                Array.Copy(mainPopHolder, innerStart, trialsPopHolder, innerStart, dim);
                            }

                        }

                        mreslis[iThread].Set();
                        //countdown.Signal();

                    } while (condition);

                    outerCountdown.Signal();

                });
            }

            countdown.Wait();

            while (popEvals < maxEvals && (Fmax - Fmin) > myEpps)
            {


                for (int i = 0; i < threadCount; i++)
                {
                    lockers[i].Set();
                }

                for (int i = 0; i < mreslis.Length; i++)
                {

                    mreslis[i].Wait();
                    mreslis[i].Reset();
                }

                double[] pom;
                pom = mainPopHolder;
                mainPopHolder = trialsPopHolder;
                trialsPopHolder = pom;

                Fmax = this.FindMax(mainPopHolder, fitnessPointers);
                Fmin = this.FindMin(mainPopHolder, fitnessPointers);

                popEvals += N;

                EvalStrategyPool(stratPtrs);

            }

            condition = false;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i].Set();
            }

            outerCountdown.Wait();

            Result result = new Result();

            result.FuncEvals = popEvals;
            result.Fmin = Fmin;
            result.Point = Common.Common.GetBestMember(mainPopHolder, fitnessPointers,N, dim);


            return result;
        }
        /// <summary>
        /// b6e6rl competitive DE, tries to find global minimum point of given function
        /// </summary>
        /// <param name="function">objective function</param>
        /// <param name="lowerBound">Lower bound vector</param>
        /// <param name="upperBound">Upper bound vector</param>
        /// <param name="maxEvals">maximum evaluations of objective function</param>
        /// <param name="dim">problem dimesion</param>
        /// <param name="N">size of population</param>
        /// <param name="myEpps">small nonnegative value, when (Fmax - Fmin) is smaller then myEpps, the termination criteria are satisfied</param>
        /// <param name="n0">controls changes in probability of strategy selection</param>
        /// <param name="delta">critial minimal value of strategy propability</param>
        /// <returns></returns>
        public Result b6e6rl(Func<double[], double> function, double[] lowerBound, double[] upperBound, int maxEvals, int dim, int N, double myEpps, double n0, double delta)
        {
            mre = new ManualResetEvent(false);
            countdown = new CountdownEvent(threadCount);
            outerCountdown = new CountdownEvent(threadCount);

            lockers = new ManualResetEventSlim[threadCount];

            mreslis = new ManualResetEventSlim[threadCount];

            condition = true;

            this.n0 = n0;
            this.delta = delta;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i] = new ManualResetEventSlim(false);
            }


            for (int i = 0; i < mreslis.Length; i++)
            {
                mreslis[i] = new ManualResetEventSlim(false);
            }

            //random stream initialization
            RngStream stream = new RngStream();
            long seed = DateTime.Now.Ticks & Int32.MaxValue;
            RngStream.setPackageSeed(new long[] { seed, seed, seed, seed, seed, seed });

            this.nRand = new Rand();

            double Fmin = double.MaxValue;
            double Fmax = double.MinValue;

            double[] uBlBDiff = new double[dim]; //upperBound - lowerBound 

            for (int i = 0; i < dim; i++)
            {
                uBlBDiff[i] = upperBound[i] - lowerBound[i];
            }

            int popEvals = N;

            int padding = 16;
            int spacing = 16;

            //pointers
            int[] populationPointers = new int[N];
            int[] fitnessPointers = new int[N];

            double[] mainPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];
            double[] trialsPopHolder = new double[padding + ((dim + 1) * N) + (threadCount * spacing)];

            int helper = 0;

            for (int threadActual = 0; threadActual < threadCount; threadActual++)
            {
                var max = N * (threadActual + 1) / threadCount;

                double[] temp = new double[dim];
                double fEvalInner = 0;

                for (int i = N * threadActual / threadCount; i < max; i++)
                {
                    //member vector initialization
                    int start = padding + i * (dim + 1) + threadActual * spacing;
                    int end = start + dim;

                    populationPointers[helper] = start;
                    fitnessPointers[helper] = end;

                    for (int j = start, forDiff = 0; j < end; j++, forDiff++)
                    {
                        temp[forDiff] = nRand.NextDouble() * uBlBDiff[forDiff] + lowerBound[forDiff];
                    }

                    fEvalInner = function.Invoke(temp); //evaluation
                    if (fEvalInner < Fmin) Fmin = fEvalInner; //is fMin ?

                    Array.Copy(temp, 0, mainPopHolder, start, dim); // temp -> mainPopHolder
                    mainPopHolder[end] = fEvalInner; // fEval -> mainPopHolder

                    helper++;
                }
            }

            Fmax = this.FindMax(mainPopHolder, fitnessPointers);
            Fmin = this.FindMin(mainPopHolder, fitnessPointers);

            // init strategy pool input variables, find real roots of polynomical equantion etc.
            // this should be a list of tuple variables F/Cr

            //get strategy F and Cr for binomical/exponantial crossover

            double[][] controlParameters = new double[12][];

            controlParameters = GetControlParameters(dim);

            //number strategy successes initialiazation/
            int longPadding = 16;
            int longSpacing = 16;

            strategySuccesses = new long[longPadding + (12 * longSpacing)];
            int[] stratPtrs = new int[12];

            for (int i = 0; i < 12; i++)
            {
                int offset = 16 + i * 16;

                strategySuccesses[offset] = 1;
                stratPtrs[i] = offset;
            }

            ResetStrategySelectionProb();
            ResetStrategySuccesses();

            // main loop
            for (int k = 0; k < threadCount; k++)
            {
                int iThread = k;

                ThreadPool.QueueUserWorkItem((st) =>
                {
                    lockers[iThread] = new ManualResetEventSlim(true);


                    int start = N * iThread / threadCount;
                    int end = N * (iThread + 1) / threadCount;
                    int n1 = N - 1;
                    int dim1 = dim - 1;
                    RngStream localStream = new RngStream();

                    double[] temp = new double[dim];
                    double innerFEval = 0;
                    int i = start;

                    //roulette
                    double rd = 0.0;
                    double sum = 0.0;

                    int strategyIndex = 0;
                    int iPtr = 0;

                    countdown.Signal();

                    do
                    {

                        lockers[iThread].Wait();
                        lockers[iThread].Reset();

                        for (i = start; i < end; i++)
                        {
                            //roulette
                            rd = localStream.NextDouble();
                            sum = 0.0;

                            strategyIndex = 0;
                            iPtr = 0;

                            do
                            {
                                sum += stratPropabilities[iPtr];
                                strategyIndex = iPtr;
                                iPtr++;

                            } while (sum < rd && iPtr < 12);

                            double F = controlParameters[strategyIndex][0];
                            double Cr = controlParameters[strategyIndex][1];

                            int r0Index = localStream.randInt(0, n1);
                            int r1Index = localStream.randInt(0, n1);
                            int r2Index = localStream.randInt(0, n1);

                            while (i == r0Index) r0Index = localStream.randInt(0, n1);

                            while (i == r1Index || r0Index == r1Index) r1Index = localStream.randInt(0, n1);

                            while (i == r2Index || r0Index == r2Index || r1Index == r2Index) r2Index = localStream.randInt(0, n1);

                            int jrand = localStream.randInt(0, dim1);

                            //RANDOM LOCALIZATION
                            int pom = 0;

                            if (mainPopHolder[fitnessPointers[r2Index]] > mainPopHolder[fitnessPointers[r1Index]])
                            {
                                if (mainPopHolder[fitnessPointers[r1Index]] < mainPopHolder[fitnessPointers[r0Index]])
                                {
                                    pom = r0Index;
                                    r0Index = r1Index;
                                    r1Index = pom;
                                }
                            }
                            else
                            {
                                if (mainPopHolder[fitnessPointers[r2Index]] < mainPopHolder[fitnessPointers[r0Index]])
                                {
                                    pom = r0Index;
                                    r0Index = r2Index;
                                    r2Index = pom;
                                }
                            }

                            int innerStart = padding + i * (dim + 1) + iThread * spacing;
                            int innerEnd = innerStart + dim;

                            int popPtr = i;

                            if (strategyIndex <= 6)//binomical crossover
                            {
                                for (int j = innerStart, jPtr = 0; j < innerEnd; j++, jPtr++)
                                {
                                    //krizeni
                                    if (localStream.randU01() <= Cr || jPtr == jrand)
                                    {
                                        temp[jPtr] = mainPopHolder[populationPointers[r0Index] + jPtr] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jPtr] -
                                            mainPopHolder[populationPointers[r2Index] + jPtr]);

                                        if (temp[jPtr] < lowerBound[jPtr] || temp[jPtr] > upperBound[jPtr])
                                        {
                                            temp[jPtr] = Common.Common.Zrdcad(temp[jPtr],
                                                lowerBound[jPtr], upperBound[jPtr]);
                                        }
                                    }
                                    else
                                    {
                                        temp[jPtr] = mainPopHolder[j];
                                    }
                                }
                            }
                            else
                            {
                                int jRandtemp = jrand;

                                if (strategyIndex > 6)
                                {
                                    for (int o = 0; o < temp.Length; o++)
                                    {
                                        temp[o] = 0;
                                    }

                                    do
                                    {
                                        // 
                                        // x = r0 + F * (r1 - r2) - mutant
                                        temp[jRandtemp] = mainPopHolder[populationPointers[r0Index] + jRandtemp] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jRandtemp] -
                                            mainPopHolder[populationPointers[r2Index] + jRandtemp]);

                                        //zrdcadleni
                                        if (temp[jRandtemp] < lowerBound[jRandtemp] || temp[jRandtemp] > upperBound[jRandtemp])
                                        {
                                            temp[jRandtemp] = Common.Common.Zrdcad(temp[jRandtemp], lowerBound[jRandtemp], upperBound[jRandtemp]);
                                        }

                                        ++jRandtemp; //increment

                                    } while (localStream.NextDouble() <= Cr && jRandtemp < dim);

                                    //pokud rand <= Cr pretece populaci
                                    if (jRandtemp == dim - 1)
                                    {
                                        jRandtemp = 0;

                                        do
                                        {
                                            // 
                                            // x = r0 + F * (r1 - r2) - mutant
                                            temp[jRandtemp] = mainPopHolder[populationPointers[r0Index] + jRandtemp] +
                                            F * (mainPopHolder[populationPointers[r1Index] + jRandtemp] -
                                            mainPopHolder[populationPointers[r2Index] + jRandtemp]);

                                            //zrdcadleni
                                            if (temp[jRandtemp] < lowerBound[jRandtemp] || temp[jRandtemp] > upperBound[jRandtemp])
                                            {
                                                temp[jRandtemp] = Common.Common.Zrdcad(temp[jRandtemp], lowerBound[jRandtemp], upperBound[jRandtemp]);
                                            }

                                            ++jRandtemp; //increment  

                                        } while (localStream.NextDouble() <= Cr && jRandtemp < jrand);

                                        if (jRandtemp < jrand)
                                        {
                                            for (int j = jRandtemp; j < jrand; j++)
                                            {
                                                temp[j] = mainPopHolder[innerStart + j];
                                            }
                                        }
                                    }
                                    else //jinak dokopirovani
                                    {
                                        for (int j = jRandtemp; j < dim; j++)
                                        {
                                            temp[j] = mainPopHolder[innerStart + j];
                                        }

                                        for (int j = 0; j < jrand; j++)
                                        {
                                            temp[j] = mainPopHolder[innerStart + j];
                                        }
                                    }//endElse
                                }//endIf
                            }//endElse

                            innerFEval = function.Invoke(temp);

                            if (innerFEval <= mainPopHolder[innerEnd])
                            {
                                Interlocked.Increment(ref strategySuccesses[stratPtrs[strategyIndex]]);

                                trialsPopHolder[innerEnd] = innerFEval;
                                Array.Copy(temp, 0, trialsPopHolder, innerStart, dim);
                            }//succesful trial
                            else
                            {
                                trialsPopHolder[innerEnd] = mainPopHolder[innerEnd];
                                Array.Copy(mainPopHolder, innerStart, trialsPopHolder, innerStart, dim);
                            }

                        }

                        mreslis[iThread].Set();
                        //countdown.Signal();

                    } while (condition);

                    outerCountdown.Signal();

                });
            }

            countdown.Wait();

            while (popEvals < maxEvals && (Fmax - Fmin) > myEpps)
            {


                for (int i = 0; i < threadCount; i++)
                {
                    lockers[i].Set();
                }

                for (int i = 0; i < mreslis.Length; i++)
                {

                    mreslis[i].Wait();
                    mreslis[i].Reset();
                }

                double[] pom;
                pom = mainPopHolder;
                mainPopHolder = trialsPopHolder;
                trialsPopHolder = pom;

                Fmax = this.FindMax(mainPopHolder, fitnessPointers);
                Fmin = this.FindMin(mainPopHolder, fitnessPointers);

                popEvals += N;

                EvalStrategyPool(stratPtrs);

            }

            condition = false;

            for (int i = 0; i < lockers.Length; i++)
            {
                lockers[i].Set();
            }

            outerCountdown.Wait();

            Result result = new Result();

            result.FuncEvals = popEvals;
            result.Fmin = Fmin;
            result.Point = Common.Common.GetBestMember(mainPopHolder, fitnessPointers, N, dim);


            return result;
        }

        private void EvalStrategyPool(int[] stratSuccsPointers)
        {
            lock (locker)
            {
                double sum1 = NjSum(stratSuccsPointers);

                for (int i = 0; i < stratPropabilities.Length; i++)
                {
                    stratPropabilities[i] = (strategySuccesses[stratSuccsPointers[i]] + n0) / sum1;

                    if (stratPropabilities[i] < delta)
                    {
                        this.ResetStrategySelectionProb();
                        this.ResetStrategySuccesses();

                        return;
                    }
                }
            }


        }

        private double NjSum(int[] stratSuccsPointers)
        {
            double sum = 0;

            for (int i = 0; i < stratSuccsPointers.Length; i++)
            {
                sum += (strategySuccesses[stratSuccsPointers[i]] + n0);
            }

            return sum;
        }

        private void ResetStrategySuccesses()
        {
            for (int i = 0; i < 12; i++)
            {
                int offset = 16 + i * 16;

                strategySuccesses[offset] = 1;
            }


        }

        private void ResetStrategySelectionProb()
        {
            for (int i = 0; i < stratPropabilities.Length; i++)
            {
                stratPropabilities[i] = 1 / H;
            }
        }

        private double[][] GetControlParameters(int dim)
        {
            double[][] controlParameters = new double[12][];

            for (int i = 0; i < controlParameters.Length; i++)
            {
                controlParameters[i] = new double[2];
            }

            double[] binF = new double[2] { 0.5, 0.8 };
            double[] binCr = new double[3] { 0, 0.5, 1.0 };

            for (int i = 0, count = 0; i < binF.Length; i++)
            {
                for (int j = 0; j < binCr.Length; j++, count++)
                {
                    controlParameters[count][0] = binF[i];
                    controlParameters[count][1] = binCr[j];
                    //pool[count].InitParams(binF[i], binCr[j], _dim, _N, _lowerBound, _upperBound, _nRand);
                }
            }

            double[] expF = new double[2] { 0.5, 0.8 };
            double[] expCr = GetExpCrValues(dim);

            for (int i = 0, count = 6; i < expF.Length; i++)
            {
                for (int j = 0; j < expCr.Length; j++, count++)
                {
                    controlParameters[count][0] = expF[i];
                    controlParameters[count][1] = expCr[j];
                    //pool[count].InitParams(expF[i], expCr[j], _dim, _N, _lowerBound, _upperBound, _nRand);
                }
            }

            return controlParameters;
        }

        private double[] GetExpCrValues(int _dim)
        {
            double[] CrVals = new double[3];

            double dim = (double)_dim;

            double p2 = (1 + 1 / dim) / 2.0;
            double p1 = (1 / dim + p2) / 2.0;
            double p3 = (p2 + 1) / 2.0;

            double[] pVals = new double[3];

            pVals[0] = p1;
            pVals[1] = p2;
            pVals[2] = p3;

            for (int i = 0; i < CrVals.Length; i++)
            {
                double[] polyCoeffs = new double[_dim + 1];

                polyCoeffs[_dim] = 1;
                polyCoeffs[0] = dim * pVals[i] - 1;
                polyCoeffs[1] = -dim * pVals[i];

                CrVals[i] = MathNet.Numerics.RootFinding.Bisection.FindRoot(x => (Math.Pow(x, dim) + polyCoeffs[0] + polyCoeffs[1] * x), 0, 0.999999,
                    1e-8, 100);
            }

            return CrVals;
        }

    }
}
